Constraint-preserving boundary conditions in the Z4 Numerical Relativity formalism 
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The constraint-preserving approach is discussed in parallel with other recent developments with 
the goal of providing consistent boundary conditions for Numerical Relativity simulations. The 
case of the first order version of the Z4 system is considered, and constraint-preserving boundary 
conditions of the Sommerfeld type are provided. The stability of the proposed boundary conditions 
is related with the choices of the ordering parameter. This relationship is explored numerically 
and some values of the ordering parameter are shown to provide stable boundary conditions in the 
absence of corners and edges. Maximally dissipative boundary conditions are also implemented. In 
this case, a wider range of values of the ordering parameter is allowed, which is shown numerically 
to provide stable boundary conditions even in the presence of corners and edges. 
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I. INTRODUCTION 

The relevance of the initial boundary value problem 
(IBVP) for Numerical Relativity has been pointed out 
many times since the ground-breaking work of Stew- 
art Q. 

The origin of the problem is the well-known fact that 
Einstein's equations 



(1) 



when interpreted as second-order field equations for the 
metric components g^ u , provide only six evolution equa- 
tions for the space components whereas the remaining 
four Einstein's equations 



(2) 



are just second-order constraints on 0]. The evolution 
equations are then a reduction of the full Einstein's sys- 
tem. Notice that this reduction is not uniquely defined, 
as far as one can add in multiples of the constraints to 
the original evolution system: this freedom is at the root 
of the diversity of the proposed evolution formalisms (for 
a review, see Ref. |j| )• 

The current approach in Numerical Relativity is to use 
one of these reductions of Einstein's field equations, plus 
four coordinate conditions, as the main evolution system 
in order to compute the full set of metric components. 
This is the unconstrained, or free evolution approach 
in which the constraint equations are mainly used 
for monitoring the accuracy of the simulations. As a 
consequence, the evolution system has an extended space 
of solutions which contains, in addition to the true ones, 
constraint- violating solutions that do not verify the full 
set of Einstein's equations. 

This poses the question of which are the requirements 
in order to get, in the free evolution approach, true Ein- 
stein's solutions in a consistent and stable way, avoiding 
any drifting towards extended, constraint-violating so- 
lutions 0. The key point is to analyze the subsidiary 



system 



V V {CT - 8tt TH = 0, 



(3) 



which follows from the contracted Bianchi identities and 
the conservation of the stress-energy tensor. As it is 
well known, the subsidiary system ensures that the con- 
straints J5J) are first integrals of the main evolution sys- 
tem. But it can also be interpreted as providing evolution 
equations for the constraints deviations. Of course, the 
same is true for any variant obtained by combining @ 
with space derivatives of either evolution or constraint 
equations: this freedom can be easily used in order to 
get a variant of the subsidiary system |J3J with a strongly 
hyperbolic even symmetric hyperbolic, principal part: 



d (G 00 - 87r T 00 ) + d k (G ok -. 



!tt T ok ) 



d (G k ° - 8tt T k °) + d k {G Q0 - 8tt T 00 ) = 



(4) 
(5) 



(normal coordinates). 

In the case of the pure initial value problem (IVP), 
the question about the consistency and stability for free- 
evolution true-Einstein's solutions has been given a pre- 
cise answer in Ref. p]j: 

• Consistency: The initial data must verify the con- 
straint equations. 

• Stability: The principal part of both the main evo- 
lution system and the subsidiary system must be 
strongly hyperbolic. 

Our former considerations suggest that the strong hyper- 
bolicity requirement on the subsidiary system is always 
fulfilled in the cases in which the only constraints are the 
energy- momentum ones @. 

In the case of the IBVP, which is the main subject of 
this paper, a number of developments are currently under 
way. Let us briefly summarize the main ones: 
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1. Constraint-preserving boundary conditions. 

The original work in Ref. [ll was focused on the 
Frittelli-Reula evolution system p> Q . More recently H, 
0,0,0], the constraint-preserving boundaries approach 
has been extended to other symmetric-hyperbolic sys- 
tems of the Kidder-Scheel-Teukolsky (KST) type [Hll3l|. 
The full programme consists in the following steps: 

• Writing down the subsidiary system (it can be © 
or any variant of it) as a first order evolution sys- 
tem for constraint deviations, with symmetric hy- 
perbolic principal part. 

• Providing algebraic (Dirichlet) boundary condi- 
tions for the incoming modes of the subsidiary sys- 
tem in such a way that the total amount of con- 
straint deviations (as measured with a suitable en- 
ergy estimate) keeps bounded. 

• Interpreting these algebraic boundary conditions 
of the subsidiary system as differential (Neumann) 
boundary conditions for the constraint-related in- 
coming modes of the main evolution system. 

• Completing the resulting subset of boundary condi- 
tions by adding suitable conditions for the remain- 
ing incoming modes. 

• Checking the stability of the final set of the main 
system boundary conditions. In Ref. [I.], the 
Majda-Osher theory is applied and the uniform 
Kreiss condition [a is obtained as a result. 



2. Einstein's boundary conditions 

An alternative approach has been proposed by Frittelli 
and Gomez [ljj [Is 13 • For a better understanding, 
let us start by writing down the constraint equations (0) 
in a covariant form 

n !/ (G^-8 7 rT^) = 0, (6) 

where n v is the normal to the constant-time hypersur- 
faces. In local adapted coordinates we have 

n v = a 8® , (7) 

so that the original form J2J) is recovered. The absence of 
second time derivatives in J2J) can now be rephrased as 
the absence of second derivatives normal to the spacelike 
constant-time hypersurfaces. 

Now let us invoke the general covariancc of Einstein's 
theory. We will realize that the same kind of result 
must hold for other kind of hypersurfaces, not just the 
spacelike ones. We can then start again from the co- 
variant form JfJJl but with n v being now the normal to 
any timelike hypersurface, like the ones corresponding to 



the boundaries of our computational domain. In local 
adapted coordinates, we could take for instance 

n» ~ SJ > (8) 

so that no second derivatives normal to the constant-z 
hypersurfaces would appear in JBJ. 

If the main evolution system is written in first order 
form, the absence of second order normal derivatives en- 
sures that four combinations of the first-order dynam- 
ical fields can be consistently computed at the bound- 
aries without any recourse to outside information. The 
Fritelli-Gomez idea is to use these combinations in order 
to get consistent boundary conditions for a subset of four 
incoming modes (Einstein's boundaries). Of course, suit- 
able conditions for the remaining incoming modes must 
also be provided and the stability of the full set must be 
checked. 



3. Harmonic coordinates 

Still another appro ach is due to Szilagyi, Winicour and 
coworkers [HJ, |2(J, 12 ll |22| ■ Their formulation looks quite 
different from the preceding ones, so that we will need 
to rephrase some statements in order to point out the 
underlying similarities. 

For instance, instead of the reduction of Einstein's 
equations, we will consider the equivalent extension of the 
solution space. In the harmonic coordinates approach, 
this extension is achieved by writing down the principal 
part of Einstein's equations as a set of generalized wave 
equations with some extra terms (de Donder-Fock de- 
composition j23l. l24l | ) and then getting rid of these extra 
terms by requiring the four spacetime coordinates to be 
harmonic functions, that is 

ax" = (9) 

(x M is considered as a set of four scalar functions here). 

Let us compare now the Harmonic Coordinates ap- 
proach with the preceding ones: 

• The resulting (relaxed) system is used as the main 
evolution system for the full set of metric compo- 
nents. By construction, its principal part amounts 
to a set of wave equations, so that symmetric hy- 
perbolicity is ensured. 

• True Einstein's solutions are recovered only when 
imposing the coordinate conditions © , which play 
here the role of the constraints. The extra prin- 
cipal terms that were suppressed from the origi- 
nal Einstein's system contained first derivatives of 
these coordinate constraints (second derivatives of 
the metric components) [2lj . 

• The subsidiary system can be again obtained from 
(J3J. Terms containing the main evolution system 
equations, or their derivatives, vanish separately so 
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that only the contribution of the extra principal 
terms remains. Notice that the resulting subsidiary 
system is of second order in the coordinate con- 
straints JJjJ, in contrast with the former approaches. 

• The principal part of the (second order) subsidiary 
system is symmetric hyperbolic: it amounts again 
to a set of wave equations on the coordinate con- 
straints ©. 

In Ref. plj . a boundary condition derived from the 
local reflection symmetry requirement is analyzed. Con- 
straint preservation is explicitly shown and a theorem of 
Secchi 25] is used in order to show that the resulting 
IBVP is well posed. This theoretical result is checked by 
means of the numerical robust-stability test (2^, which 
is adapted so that reflection boundary conditions are ap- 
plied along just one space axis, while keeping periodic 
boundary conditions along the other two. 

Due to the limited use of reflection symmetry bound- 
ary conditions in practical applications, a proposal is 
made for extending these results to boundary condi- 
tions of the Sommerfeld type, as it was done previously 
in a different framework 1271 . 



4- The Z4 case 

The Z4 approach pHI^ uses an extra dynamical four- 
vector, along the track of previous formulations which 
contained extra dynamical quantities [23, [^J [53, US • It 
has strong similarities with the harmonic coordinates ap- 
proach, while providing much more flexibility regarding 
the gauge choices. 

• The main evolution system is obtained by modify- 
ing Einstein's equations with the help of the extra 
four- vector Z^, namely 

GW + + v - z m - (V P Z p ) g„ v = 8nT^. (10) 

This system provides ten evolution equations for 
the set formed by the six space components of the 
metric plus the four Z^ components. A first or- 
der version can be easily obtained which, when 
supplemented with suitable gauge conditions, has 
been shown to have a strongly hyperbolic principal 
part [H. 

• True Einstein's solutions can be recovered by re- 
quiring the vanishing of the extra four-vector 

Z" = 0, (11) 

so that this condition can be considered as a set 
of four algebraic constraints. Notice that the main 
evolution system (jlUf) is of a mixed type: it contains 
second order derivatives of the metric, but only first 
order derivatives of the extra four-vector. 



• The subsidiary system can be obtained from the 
covariant divergence of the main system (TQ) : h 
will be of third order in the metric and of second 
order on the constraint variables Z^. Allowing for 
(j2J), the Einstein's tensor contribution vanishes sep- 
arately, so that only the contribution of the extra 
terms remains, namely 



V„ [ VZ U + VZ" - (V P ZP) ] = , (12) 



which can be also expressed in the equivalent 
form 



□ Z^ + R^Z V = . (13) 



• It follows from (|13fl that the subsidiary system, 
when considered as a second order system for the 
algebraic constraints deviations Z^ , has a symmet- 
ric hyperbolic principal part: it amounts here again 
to an uncoupled set of wave equations. 



The fact that the subsidiary system is of second or- 
der means that the vanishing of both Z^ and its first 
time derivative must be imposed on the initial data if 
one wants to ensure a priori that the resulting solution 
will be a true Einstein's one. This amounts to impose 
the usual energy and momentum constraints © on the 
initial data hypersurface, so that it can seem that the Z4 
formalism is not being of much help. But on the other 
side, if one is checking a given solution a posteriori, the 
vanishing of Z\i in a given spacetime domain ensures that 
the same is true for its derivatives, so that this solution 
is necessarily a true Einstein's one. This is why one can 
monitor constraint violations by looking just at the val- 
ues of Z^ and, more important, this is why one can devise 
constraint-preserving strategies by aiming at the vanish- 
ing of Z^. Here is where the Z4 formalism shows its main 
advantages. 

In what follows, we will consider the fully first-order 
version of the Z4 system |29f , as summarized in Section 
2. In Section 3, we will apply the constraint-preserving 
boundary conditions programme, obtaining as a result 
conditions of the Sommerfeld type for the main evolu- 
tion system. The stability of these conditions is studied 
in Section 4, including the use of the robust-stability nu- 
merical test. As a result, our Sommerfeld-like conditions 
will be shown to behave in the same way as the reflection 
symmetry ones proposed in Refs. [2lll22]|. 
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II. FIRST ORDER Z4 SYSTEM 

The general-covariant equations (|10fl can be written in 
the equivalent 3+1 form [28| | 



(d t - Cp) 
(d t - Cp) K i:j 



(d t - Cp) © 
(d t - Cp) z t 



-2 a Kj 



-ViQij 



(14) 



[ {3) R, 



2Kl + (trif-20) K K 
Sij + - (tr S - t) jij } 



(15) 



| [ (3) i? + 2 V k Z k + (trK-2e)trK 

- tr(K 2 ) - 2r ] - Z k a k (16) 
a [ Vj (i^ - V'trif) + <%© 

(17) 



2 -Ki^ Z q — Si 



e. 



where we have noted 



r = 8-ko 2 T 00 , ^ = 8ira T% , ee 8vr T, 



(18) 



In the form I|14ll7fl . it is evident that the Z4 evolu- 
tion system is fully relaxed: it consists only of evolution 
equations. The original constraints 1|11|> . which can be 
translated into 







0. 



Z, 



0, 



(19) 



are algebraic so that the full set of field equations 11U|) 
is actually used during evolution, like in the harmonic 
coordinates case. 

But now we have not to impose the harmonic coor- 
dinate conditions ©. We will consider instead a wider 
class of gauge conditions, in which the time slicing will 
be of the form (23 



(d t - Cp) In a = - fa (tr K - m0) 



(20) 



(generalized harmonic sli cing ). Although more general 
cases can be considered |34j . we will use here normal 
coordinates (zero shift) for simplicity. 

A first order version of the Z4 evolution system (|14I17|) 
can be obtained by introducing the first space derivatives 



Ak = atk/a, D ki: j = - d^ij (21) 

as independent dynamical quantities, so that the full set 
of dynamical fields can be given by 

u = {a, 7l y, Kij, A k , D Hj , 0, Z k } (22) 

(38 independent fields). 

Of course, one must provide evolution equations for 
the new quantities 1)21(1 : the simplest way is to take 

d t A k + d k [fa(trK-m&)] = (23) 
d t D kij + dkiaKij ] = . (24) 



Notice that one could add to 1(231 I24fl a number of terms 
involving first derivatives of either or Z k . This would 
amount to introduce coupling terms with either the En- 
ergy or the Momentum constraints, as in the KST sys- 
tem |l2t 1 1 3| . each one with its own free parameter. 

We have chosen instead to keep the simplest form ()23l 
I24JI because the first order constraints 1)21(1 evolve in a 
trivial way, that is 



d t [ A k - d k In a } = 
d t [ D kij - d k 7y ] = 



(25) 
(26) 



so that the relationship between the first and the second 
order versions of the evolution system is more transpar- 
ent. We are losing in this way the possibility of playing 
with a number of extra free parameters. 

Care must be taken, however, when expressing the 
Ricci tensor ^'Rij in (|15|l in terms of the derivatives of 
Dkij , because as far as the definitions 1(21(1 are no longer 
enforced, the identities 



C kl = d [k A l} = 



Ckuj = d[ k Dqij — 



(27) 



can not be taken for granted in first order systems. As 
a consequence of these ordering ambiguities, the princi- 
pal part of the evolution equation 1(15(1 leads to a one- 
parameter family of non-equivalent first-order versions, 
namely 



dtKn + d k [ a \% 



(28) 



where 



D 



Dji - SfE 3 - 5 k E t 



1 ak/ A 7~i i it/ ^ i 1 ski 



+ - Sf(A, - Dj + 2V j ) + - 5](Ai - A + 2Vi) (29) 

and we have noted 

Di = Y s D zrs , Ei ee y rs D rsi , V k ee D k -Ek-Z k . (30) 

Notice that the parameter choice £ = +1 corresponds to 
the standard Ricci decomposition 

^ = dk r ij &i r kj ~\~ r rfeT y r r iV kj (3i) 

whereas the opposite choice £ = — 1 corresponds to the 
de Donder-Fock [23ll24| decomposition 



d k D k v] + d (l T. l)k k - 2D/ k D 



4-D iD rs j T^ij-gJ^j 



kij 



r r r 



(32) 



which is most commonly used in Numerical Relativity 
formalisms. The ordering ambiguities do not affect to 
the principal part of eq. ((161) . namely 

d t e + d k [aV k ] = ... (33) 

The resulting first order system has been shown to be 
strongly hyperbolic (2i| provided that the first gauge pa- 
rameter / is greater than zero. In the harmonic slic- 
ing case (/ = 1), the second gauge parameter is fixed 
(to = 2). The full list of eigenvectors is given in Ap- 
pendix A. 
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III. CONSTRAINT-PRESERVING BOUNDARY 
CONDITIONS 

We have seen in the Introduction that the simple equa- 
tion H13|) provides the subsidiary system for the devia- 
tions of the algebraic constraints This would be 
the whole story if we were planning to use the second 
order version i|l(J|) of the evolution system. But we prefer 
to focus here in the first-order-in-space version, as de- 
scribed in the previous section. The reason is that the 
mathematical theory of first order systems seems to be 
more developed, both at the continuum and at the dis- 
crete level, so more powerful tools are available: Energy 
methods, Total- Variation-Diminishing algorithms and so 
on (see for instance Refs. |al35|). 

There is a price to pay for this. We have found in the 
previous Section new constraints, like (|T7Jl . arising from 
ordering ambiguities in the space derivatives. The or- 
dering parameter £ appeared precisely from the coupling 
of these ordering constraints with the evolution system. 
The original subsidiary system 1(13(1 must then be ex- 
tended in order to include both these coupling terms and 
the evolution of the ordering constraints themselves. 



1. First- order subsidiary system. 

The easiest way of obtaining the full subsidiary sys- 
tem in the first-order case is just by computing the time 
derivative of the full Z4 first-order system. We give here 
(the principal part of) the resulting subsidiary system 

d t C kl = (34) 

d t C klij = (35) 

l/a 2 d 2 t 6 - A 6 = • • • (36) 

l/a 2 d 2 t Zt -AZi = j kl d k [ C a + 7 rs ( Cars 

+ (C - 1) C rUi + (C + 1) C risl )] + ••• (37) 

(the dots stand for non-principal terms). 

The subsidiary system (I34l - I37|l can be put in first order 
form in the usual way, by considering the first derivatives 
of (0, Zi) as new independent variables. The following 
evolution conditions 

d t (d k &) -d k [d t e] = o (38) 

d t (d k Zi) - d k [ d t Z t ] = (39) 

could be added then to complete (the first order version 
of) the subsidiary system. 

Notice that the evolution equations (|34l I35fl for the 
ordering constraints are trivial. This means that the 
ordering constraints themselves are eigenfields of the 
full subsidiary system (|34l - I39fl with zero characteris- 
tic speed. Moreover, the evolution equations (|36l I38|) for 
(the derivatives of) form a separate subsystem with the 
structure of the wave equation. Concerning the remain- 
ing equations (|37l I39fl . one can express them in terms of 



the quantities 

Z k i = d k Zi + dk + 7 rs [ Ci krs 

+ (C - 1) C rksi + (C + 1) C nsk ] , (40) 

so that they read 

l/a 2 d t (dtZi) ~ d k [ Z kt ] = • • • (41) 
d t Z kl -d k [d t Zi] = ■■■ , (42) 

and we get again the structure of the wave equation. 

It follows that the principal part of the subsidiary sys- 
tem (|34l - I39fl can be put in symmetric hyperbolic form. 
The characteristic speeds are either zero or the light 
speed. A simple energy estimate is provided by 

e = i/a 2 [ (d t e) 2 + ^ (d t z t )(d t z 3 )} 

+ (d k O)(d k Q) + ZyZ*. (43) 

We are now in position to take the second step in the 
constraint-preserving boundary conditions programme. 
We will impose the vanishing of all the incoming modes 
of (|34I - I39|) at the boundaries, that is 

l/a d t + n k d k 9 = (44) 
l/a d t Z t + n k Z kl = , (45) 

where ft stands here for the outwards-pointing unit nor- 
mal to the boundary surface. 

Equations 1(441 145(1 meet the two requirements we were 
looking for: 

• They provide maximally-dissipative algebraic 
boundary conditions for the subsidiary system 1(341 
I39() . In this way, no constraint-violating modes 
are allowed to enter across the selected boundary. 

• They will provide, as we will see in what follows, 
four boundary conditions of the Sommerfeld type 
for the evolution system 1(14117(1 . which can be con- 
sistently imposed in order to obtain true solutions 
of Einstein's field equations. Notice that the extra 
terms in the definition ((40(1 consist in ordering con- 
straints, which would not appear in a second order 
in space formulation. 



2. Boundary conditions implementation. 

The third step in the programme is to use the resulting 
values (6 (&otm) , zf oun) ), as computed from gUHUl, in 
order to obtain four of the main system's incoming fields 
at the boundary. This process is not free from ambigui- 
ties, like the choice of a suitable basis for the dynamical 
fields. 

For a symmetric hyperbolic evolution system, one 
could find a (positive definite) quadratic form which 
would provide a metric for the space of dynamical fields. 
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The natural choice would be then to build an orthogo- 
nal basis of dynamical fields containing both and Zi 
(or some equivalent combinations). Imposing boundary 
conditions would then consist in prescribing the values 
(|44l I45|) for these fields, while leaving the remaining ones 
unchanged. 

But the evolution system (|14I17(I is not symmetric hy- 
perbolic. This means that we do not have a unique pre- 
scription for imposing the boundary conditions, as far 
as we have many ways of selecting an appropriate set of 
dynamical fields at the boundary. A convenient starting 
point in this case is to replace the original basis 

(9, K i3 , Zi, D kij , A % ) (46) 

by one which is more adapted to the characteristic de- 
composition at the boundary, namely 

(6 , Kij , Z, , D Uj , D n±± , Vi , A , A) , (47) 

where the symbol _L replacing an index means the pro- 
jection orthogonal to ft. We have noted as D n ±± the 
traceless part of D n ±± and 

Kij = K l0 - - 7y . (48) 

Notice that the quantities D nn ±, tr(D n ±±) do not ap- 
pear explicitly in the new basis. These components must 
be computed instead from (Z + V)i and the D±ij com- 
ponents. Allowing for the definition l|H()|l. we actually 
get 

D nn± = D ± -h rs D rs± -(Z + V)± (49) 
h rs D nrs = h rs D rsn + (Z + V)» , (50) 

where h rs stands for the (inverse) metric on the boundary 
surface, namely 

h rs = 7 rs - n r n s . (51) 

The new basis (14 711 has been chosen in such a way that, 
as we can easily verify, the values of (8, Zi) appear in 
only eight eigcnficlds (four characteristic cones), namely: 

E ± = Q±V n , (52) 

L± ± = K n± ± [ X - (A ± + D ± - 2 Z ± ) 

- ^ D ±nn + ^ h rs D rs± } , (53) 

L± = h rs K rs ± [ Z n - C h rs D rsn } , 

where we have noted 

L ± = h rs L± - E ± . (54) 

In order to set up the required four boundary con- 
ditions, we will simply replace the original values for 
(6 ,Z-) by (e {bound ^ ,zf ound) ), while leaving the other 
fields in the basis (|47|l unchanged. To be more specific: 



• The original values for (9 , Zi) are replaced by 

( (boW) ^ z (bound)^ ag computed from g^gT^ re _ 

spectively. This amounts, modulo some linear com- 
binations with tangent fields (transverse deriva- 
tives), to prescribe the first term in and the 
second terms in (£ nJ _, L^). 

• The values of their 'counterpart' fields 
{V n , K n ± , h rs K rs ) are not changed by the 
boundary conditions. 

It is clear then that the values of the four characteristic 
cones (|52l 1531 154fl have been prescribed in such a way 
that the four equations 1441 145|) hold true at the selected 
boundary. 

A further source of ambiguity comes from the prescrip- 
tion of the remaining incoming eigenfields (the gauge and 
the transverse traceless ones). We will use here a conve- 
nient generalization of the maximally dissipative bound- 
ary conditions, namely 

d t G-=Q, d t [Ll±-\ (h rs L- s ) 7±± ] = 0, (55) 

although we are aware that more sophisticated choices 
could be required in physical applications. 



IV. CONSTRAINTS STABILITY 

The final step in the proposed programme is to check 
the stability of the constraint-preserving boundary con- 
ditions P H 1151155)1 . 

Notice however that the main evolution system is just 
strongly hyperbolic, but not symmetric hyperbolic (at 
least not in the generic case |28|). This means that the 
Majda-Osher theory 0] can not be directly applied, and 
the same is true for the Secchi theorems |25|. This is 
why we will check the stability of l|44l 1451 155(1 by other 
methods, both at the theoretical and the numerical level. 

From the theoretical point of view, the well-known 
Fourier-Laplace method y| could provide necessary con- 
ditions for stability ^(j- We will prefer here a simpler 
approach, by analyzing the system of equations verified 
by the dynamical fields at the boundary. We will call it 
the modified system in order to distinguish it from the 
original evolution system, which is being used at the inte- 
rior points. We will see that this approach provides some 
insight about the behavior at the boundary points. The 
drawback is that boundary points form just the outer- 
most layer of the computational domain. It follows that 
the modified system analysis has to be considered at this 
stage just as an heuristic approach, so that the stability 
of the boundary conditions must be confirmed by other 
means. More details are provided in Ref. 36]. 
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A. The modified system approach 

For the sake of clarity, let us focus first on the subset 
of dynamical fields spanned by (0, Vi). As stated in 
the previous Section, the boundary conditions are not 
affecting any of the Vi components. This means that the 
boundary values of V, verify the main evolution system 
equations, namely 



i/a d t v + d t e = 



(56) 



The original equation (|lufl for 9, however, no longer 
holds at the boundary, where one is imposing instead the 
advection equation (|44() . This means that, even at the 
continuum level, the evolution system is being modified 
at the boundaries. The modified system for the subset of 
dynamical fields (0, Zi) is given by 1)441 1561) . 

The modified subsystem (|44l I56f) has real non- negative 
characteristic speeds along any direction r, oblique to ft. 
They are actually 



{0, a (n ■ r) } . 



(57) 



It follows that (the principal part of) the modified sub- 
system (|44l 15611 can be interpreted on physical grounds 
as describing the outwards propagation of both and V, 
at the boundary. 

We can push one step further our analysis by consider- 
ing the particular case in which r is tangent to the bound- 
ary, that is orthogonal to ft. In this case the speeds (|57|) 
are fully degenerate, and a non-diagonal coupling term 
remains in (|56|l . so that the modified subsystem is just 
weakly hyperbolic. This has some relevant consequences. 
Let us assume for instance that ft is aligned with the x 
coordinate axis and that we get a static profile for of 
the form 



@ = g(y, z) 



(58) 



which trivially satisfies equation (|44|l . The derivative 
coupling in (|56|l allows then modes in the V y and V z com- 
ponents which grow in time in a linear way. These lin- 
early growing modes will actually show up in numerical 
tests, as we will see below. 

The analysis of the full modified system can be simpli- 
fied by writing down (the principal part of) the modified 
evolution equations for the combinations corresponding 
to incoming modes of the original system. Allowing for 
1511. we have 



(59) 
h rs d r [ D ns± - D sn± + ^=-1 K.±] 

d± [ K nn - f -^~ trK + A n 

+ K- (3 - 2 7 )/ + 1 0] (60) 

h r [ Dnns ~ -&snn C ^-ns 

+ A s + V s ] , (61) 



1/a d t E~ 
1/a d t L~ s 



1/a d t LT 



plus the trivial evolution equations l|55|) for the gauge and 
the transverse traceless incoming modes. 

Notice that only derivatives tangent to the boundary 
appear on the modified system equations 1)591- 1611) for the 
incoming modes. This means that all the characteristic 
speeds along the longitudinal direction ft are real and 
non-negative: they are actually 



0, a, aV?- 



(62) 



The corresponding eigenvectors are either standing fields 

( v = oy. 

A±, D ±ij , Ak-fDk + fmVk, E~ , Lr. , G~ (63) 
or outgoing fields (v — a, a*JJ): 



j Zi , Zfl 



G 



(64) 



These fields span the whole dynamical space; the modi- 
fied system is then strongly hyperbolic along the direction 
ft normal to the boundary. 

Computing the characteristic speeds along a generic 
direction r, oblique to ft, and for an arbitrary value of the 
ordering parameter, is a much harder task, even using an 
algebraic computing program. We have just checked the 
particular cases 



C = 0, ±1 



(65) 



and we have found that the modified system is at least 
weakly hyperbolic (real characteristic speeds) only in the 
C = case. This suggests that the ( = case, corre- 
sponding to a symmetric ordering of the space deriva- 
tives, could be free of boundary instabilities, as we will 
confirm below. 



B. The robust stability test 

The robust stability numerical test |2^| amounts to 
consider small perturbations of Minkowski space-time 
which are generated by taking random initial data for 
every dynamical field in the system. The level of the ran- 
dom noise must be small enough to make sure that we 
will keep in the linear regime even for hundreds of cross- 
ing times (the time that a light ray will take to cross the 
longest way along the numerical domain). We are tak- 
ing advantage in this way of the peculiar nature of the 
Einstein's equations, where the principal part is quasilin- 
ear and the non-principal (source) terms are quadratic in 
the dynamical fields. Checking the linear regime of Ein- 
stein's equations amounts then to test the behavior of 
their principal part. 

This test has been previously used [29] for to check 
numerically the stability properties of the Z4 evolution 
system for interior points. In order to avoid boundary 
effects, the grid had the topology of a three-torus, with 
periodic boundaries along every axis. We will now open 
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0.1 1 10 100 0.1 1 10 



crossing times crossing times 



FIG. 1: Loo norm of trK for different values of the ordering 
parameter C, and constraint-preserving boundary conditions 
along one single direction (periodic boundaries along the other 
two). The values of C, are shown with an interval of 0.5 for 
the sake of clarity, although the survey has been made with 
a finer interval of 0.1. The stable range is given by £ values 
in the interval [—0.5,0]. The decreasing of the norm in the 
stable regions is due to the dissipative boundary conditions 
1551 1 for the gauge modes 



the x faces and impose the constraint-preserving bound- 
ary conditions there, while keeping periodic boundary 
conditions along the other two axes. 

We show in Fig. ^ the L M norm of trK for different 
values of the ordering parameter Q. A spacing A£ = 0.5 
is used in the plot for the sake of clarity, although the 
numerical survey has been made with a finer spacing 
of A£ = 0.1. Our results show that the constraint- 
preserving boundary conditions 1441 145fl are stable if and 
only if £ is in the range [—0.5,0]. The behavior is the 
same in all the stable regions: the different values of £ 
just determine when the instabilities (if any) will appear. 

Notice that the robust stability analysis predicted the 
arising of linearly growing modes related to the trans- 
verse Vi components. In terms of the original basis, one 
can expect to see these modes in the quantities D xxy , 
D XX z, which are derived from V x , V z by the relationships 
(|49l I50fl . respectively. We can see for instance in Fig. [21 a 
growing linear mode in the norm of D xxy . This con- 
firms that the modified system analysis can be useful to 
anticipate the behavior of the boundary conditions under 
the robust stability test. Further evidence in this direc- 
tion is provided in the Appendix B, where maximally 
dissipative boundary conditions are considered. 

The robust stability test is also useful for checking the 
constraint-preserving character of the proposed bound- 
ary conditions 1441 14511 . As far as true Einstein's solu- 
tions are actually recovered by setting both <d and Zi 
to zero, the values of these quantities can be considered 
to be good indicators of constraint violations. We can 
monitor the norm of these quantities to check whether 
constraint violations are being injected into the compu- 



FIG. 2: Same as Fig. Q but now for the norm of D xxy 
and ( — 0. We plot here the results corresponding to three 
different resolutions, focusing on the first 10 crossing times. 
Notice that this is a logarithmic plot, so that the resolution- 
independent linear growing with unit slope corresponds actu- 
ally to a linearly growing dynamical mode. 



tational domain through the open boundaries. 

We can see in Fig.GHthat this is actually not the case. 
The values of and Zf are not growing at all, contrary to 
what happens to the D xxy components, as seen in Fig. [21 
Moreover, their norm is diminishing: we can understand 
this decreasing by noticing that the boundary conditions 
(|44l I45|) are, modulo some coupling with ordering con- 
straints, advection equations. This means that the val- 
ues of (|44l 145(1 are just flowing out of the computational 
domain through the open boundary. 
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10 



in 



10 



0.1 



1 10 

crossing times 



FIG. 3: Same as in the previous figures, but now the norms of 
<d and Z x are plotted in order to monitor constraint violations. 
No growth can be seen, confirming that constraint violations 
are not being injected through the boundary. 
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V. DISCUSSION AND OUTLOOK 

One can wonder about wether the multi-dimensional 
character of the problem is lost when applying boundary 
conditions to just one face. This is not the case: the 
x = constant boundary surfaces are actually surfaces, 
not just points, so that oblique modes are still present 
and can lead to instabilities, as it actually occurs when 
Q = ±1. The only thing we are avoiding in this way 
is to apply constraint-preserving boundary conditions to 
corner points (which are assigned here to the y and z 
boundary surfaces, where periodic boundary conditions 
are applied). If we try to apply conditions ^441-145(1 along 
every space direction, then instabilities appear even for 
the C = choice of the ordering parameter. 

One could argue as well that the difficulties with cor- 
ners and edges could be related with the fact that the 
main evolution system has just a strongly hyperbolic, 
but not symmetric hyperbolic, principal part. This is not 
the case, as it is shown in Appendix B, where boundary 
conditions of the maximally dissipative type are success- 
fully implemented and applied to all the boundary points 
of a cartesian-like numerical grid, including corners and 
edges. 

Our results are very similar to those of Ref . [2l| , where 
the same test is applied to reflection boundary condition: 
a stable linear-growing mode is detected, which becomes 
unstable only when the boundary conditions are applied 
also to the other faces, so that the numerical grid gets 
corners and edges. Our work can be then understood as 
an extension (in a different formalism) of the results of 
Ref. [2ll to boundary conditions of the Sommerfeld type. 

In our opinion, the main problem at corner points 
comes from the inconsistency inherent to the choice of 
a (unique) normal direction there. Different faces get 
different normal vectors, but corner points belong to two 
different faces at the same time. This is not just a theo- 
retical caveat: corners and edges pose a real problem in 
practical applications, where more work should be done 
along any of the following lines: 

• Devising an specific treatment for corner points. 
The correct implementation of constraint- 
preserving boundary conditions in the presence of 
corners is still a unsolved issue. For symmetric 
hyperbolic evolution systems, using finite differ- 
ence operators satisfying the summation by parts 
rule with respect to a diagonal scalar product 
leads to stable schemes with maximally dissipative 
boundary conditions [371II1]- In our case, with a 
strongly hyperbolic evolution system, we have got 
the same result, as presented in Appendix B. 

But these results do not extend to constraint- 
preserving boundary conditions. A major difficulty 
is that compatibility conditions between boundary 
data at adjacent faces need to be satisfied if one 
wishes to obtain smooth solutions. Necessary con- 
ditions for continuous solutions have been derived 



in Ref. [9( for a symmetric hyperbolic evolution 
system. But more conditions are needed in order 
to obtain smooth solutions. Compatibility issues 
are also present at corner points between initial 
and boundary data (see for instance Chapter 9 in 
Ref. H|). 

• Building numerical grids with smooth boundaries 
(without corners and edges), so that the constraint- 
preserving boundary conditions (|4*4"l - l4"5|) can be ap- 
plied consistently in an stable way, as we have con- 
firmed numerically by means of the robust stability 
test. The construction of 'Multi-patch' numerical 
schemes, which would allow for smooth boundaries, 
has become a major research to pic in Numerical 
Relativity. See for instance Refs. |39L l40| . 



Appendix A: Characteristic decomposition of the Z4 
first order system 

Let us consider the propagation of perturbations with 
wavefront surfaces given by the unit (normal) vector n^, 
we can write the principal part of the Z4 first order sys- 
tem in matrix form 

- dt u + M n k d k u= ... , (A.l) 

a 

where u is the full array of dynamical fields (1221 . Notice 
that derivatives tangent to the wavefront surface play no 
role here. 

A straightforward analysis of the characteristic matrix 
M provides the following list of eigenfields [29| : 

• Standing eigenfields (zero eigenvalues) 

a, 7 y, A±, D ±ij , A k - fD k + fmV k , (A.2) 

where the symbol _L replacing an index means the 
projection orthogonal to 

D±ij ■. = D kij - n k n r D rij . (A. 3) 

• Light-cone eigenfields (eigenvalues ±a) 

L ± ij = [Kij — mrij tr K ] 

± [A™y - mrij tr A"] (A.4) 
E ± ee 9 ± V n , (A.5) 

where the symbol n replacing the index means the 
contraction with nt 

Kj = n k \% V n ee n k V k . (A.6) 

• Gauge eigenfields (eigenvalues iav7) 

G ± = v7 [ tr K - n 9 ] ± [ A n + (2 - fx) V ] (A.7) 
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where we have noted for short 
fm-2 



/-I 



(A.8) 



In the degenerate case (f=l), one must have m = 2, 
so that the value of \x is not fixed. The degeneracy 
allows for any combination with (|A.5|1 , as expected. 



Appendix B: Maximally dissipative boundary 
conditions 

A convenient generalization of the maximally dissipa- 
tive boundary conditions can be implemented by just im- 
posing the vanishing of (the time derivatives of) all the 
incoming modes, that is 



dt E~ = dt L- 



d t G~ 







(B.l) 



The principal part of the modified system is the much 
simpler that in the constraint-preserving case. It is, by 
construction, strongly hyperbolic along the direction ft 
normal to the boundary, with characteristic speeds given 
again by JB^J. 

We will compute the characteristic speeds along a 
generic direction r, oblique to ft, where the vector r is 
related with n by 



n cosp + s sirup , 



and we have taken 
n 2 = 



,s 2 = l. 



ft • s = . 



(B.2) 
(B.3) 



The hyperbolicity requirement amounts to demand that 
all the resulting characteristic speeds be real for any value 
of the angle ip. 

The trivial equations (|B.1|) provide 7 (remember that 



is traceless) standing eigenfields (zero characteristic 
speed) of the modified system. Another set of 17 standing 
eigenfields is given by 



D 



py > 



A k - fD k + fmV k 



(B.4) 



where p is the direction orthogonal to both vectors n 
and s. 

The remaining 14 dynamical fields can be grouped into 
the following sectors: 

• Energy sector {E + , V s }. The corresponding evo- 
lution equations are (principal part only): 



- d t V s 

a 



- d t E+ 



-sirup d s O = — — sirup d r E^ 

-d r [ V r + O cosip ] 

-d r [ E + cosip + V s sirup + • ■ ■ 



where the dots stand for coupling terms with the 
standing eigenfields (which are irrelevant for the 



eigenvalues calculation). It follows that the char- 
acteristic speeds are given by the solutions of the 
algebraic equation 

A(A — a cosip) — — a 2 sin 2 tp , (B.5) 

so that real characteristic speeds are obtained for 
every value of tp. 

• Gauge sector {G + , A s }. The corresponding evo- 
lution equations are (principal part only): 

— dtA s = —sirup d s \ fCtrK — m 0)] 
a 

= ~ 2 V7 sin( P dr G + + 

- d t G+ = -d r \ J]A r + f cosp trK ] 
a 

= — \fj d r [ G + cosip + A s sirup + ■ • • ] 

where the dots stand for coupling terms with the 
previous sectors. It follows that the characteristic 
speeds are given by the solutions of the algebraic 
equation 

A(A — a y/f cosip) — -fa 2 sin 2 ip , (B.6) 

so that, allowing for the positivity of the gauge pa- 
rameter /, real characteristic speeds are obtained 
again for every value of ip. 

• Metric sector {L^ , D s ij}. The corresponding 
evolution equations can be written as (principal 
part only) 

1 

2 



1 








dt 


D s ij 


a 




1 

o 


dt 


hi 



-sirup d s Kij — — — sirup d r [ 



(B.7) 



-d r [ Xij + cosip Ki 3 + ■■■ 

{r%Kjij ~t~ TjKfii riiK r j rijK r i) J 

-d r [ Ljf cosip + Dsij sirup H 

simp {siK n:j + SjK ni - niK SJ - n K S i) 
sirup (D ijs + Djis - SiEj - SjEi) ] . 



1 + C 
2 

i + C 

2 

i + C 



where the dots stand again for coupling terms with 
the previous sectors. 



The evolution equation i|B.8|) for these outgoing 
'metric' fields contains (unless £ = — 1) crossed coupling 
terms that complicate the analysis. One gets three vari- 
ants of the same algebraic equation 



A(A — a cosip) — — a 2 sin 2 ip 



(B. 



A(A - a cosip) = l a 2 sin 2 ip [1 - (^-^) 2 ] (B.9) 



A(A - a cosip) = 2 V t 1 ~ ^ + 01 > ( B - 10 ) 
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FIG. 4: Same as Fig. □ but now with maximally dissipative 
boundary conditions enforced along every axis. Again, the 
values of ( are shown with an interval of 0.5 for the sake 
of clarity, although the survey has been made with a finer 
interval of 0.1. It confirms that the £ < choices are stable, 
as expected from the modified system analysis. Notice that, 
in the stable cases, the decrease after 100 crossing times is 
more than one order of magnitude greater than in Fig. Q 
This shows the effect of the maximally dissipative boundary 
conditions: they are actually dissipating all the dynamical 
fields. 



depending on the particular set of components considered 
(the last two equations appear twice, so that one gets 
10 characteristic speeds that complete the full set of 38). 
The most restrictive is the last one (|B.10ll : it implies that 
we get complex characteristic speeds for some values of 
ip unless 



C<o, 



(B.ll) 



so that the standard ordering case (£ — + 1) is excluded. 

We can check out these results by using again the 
robust stability test-bed. We will enforce the maxi- 
mally dissipative conditions (|B.1|) along every axis in 
a cartesian-like numerical grid, including corners and 
edges. We will survey the values of the £ parameter in 
the interval [—1, 1], with a spacing A( = 0.1. 

We show in Fig. 21 the time evolution of the maxi- 
mum of the absolute value of trK (a spacing AC = 0.5 
is used in the plot for the sake of clarity). Our results 
show that the positive choices of the ordering parameter 
are actually unstable, whereas the choices in the range 
[—1,0] are stable and behave in the same way. Notice 



that the norm of trK in Fig. 0] is decreasing, in the sta- 
ble cases, much faster than in the constraint-preserving 
case (Fig. [TJ. This can not be explained just by the fact 
that boundary conditions are now being applied along 
the three coordinate axes: the boundary conditions are 
actually dissipating all the dynamical fields. 

We show in Fig. [3] the time evolution of the maximum 
of the absolute value of both and Zi for the symmetric 
ordering (£ = 0) case. As far as their values are di- 
minishing, one can conclude that no constraint- violating 
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FIG. 5: Same as in the previous figure, but now the norms 
of and Z x are plotted in order to monitor constraint vi- 
olations. Their decay indicates that the constraint-violating 
modes are diminishing, even much faster than in Fig. [3] But 
now the reason is a different one: the boundary conditions 
are dissipating all the dynamical fields. 



modes are being produced at the boundaries. But notice 
that this is at the price of dissipating all the dynamical 
fields. One can not conclude then that maximally dis- 
sipative boundary conditions are constraint-preserving: 
constraint-related fields are flowing out through the 
boundaries in the same way as the remaining degrees 
of freedom. 
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